A preliminary, observational study using whole-blood RNA sequencing reveals differential expression of inflammatory and bone markers post-implantation of percutaneous osseointegrated prostheses

Aims While the benefits of direct skeletal attachment of artificial limbs are well recognized, device failure due to infection and insufficient osseointegration remain obstacles to obtaining consistently successful outcomes. Currently, the potential for device failure is assessed by subjective pain, clinical function scores, radiographic evidence of bone atrophy, and the presence of radiolucent lines at the bone-implant interface, and subjective pain and function scores. Our hypothesis is that measurable biological indices might add another objective means to assess trends toward bone and stomal healing. This longitudinal cohort study was undertaken to identify potential serological biomarkers suggestive of bone remodeling and the presence of stomal tissue inflammation. Methods Ten unilateral transfemoral amputee veterans, who were implanted with a percutaneous osseointegrated (OI) skeletal limb docking system, were recruited to participate in this IRB-approved study. Venous blood samples were obtained from before the Stage 1 Surgery up to 1 year following the Stage 2 Surgery. Whole-blood RNA was extracted, sequenced, mapped, and analyzed. Of the significant differentially expressed (DEGs) genes (p<0.05) identified, four genes of interest (IL12B, IL33, COL2A1, and SOST) were validated using qPCR. Enrichment analysis was performed to identify significant (p<0.01) Gene Ontology (GO) terms. Results Most differentially expressed genes were only detected at PoS1 immediately after the first surgery. Of the significant genes identified, IL12B and IL33 were related to inflammation, and COL2A1 and SOST were associated with bone remodeling. These four genes were identified with greater than 20 log fold-change. Conclusion Whole-blood RNA-seq data from 10 patients who previously underwent percutaneous osseointegrated lower limb implantation revealed four genes of interest that are known to be involved in inflammation or bone remodeling. If verified in future studies, these genes may serve as markers for predicting optimal bone remodeling and stomal tissue healing following OI device implantation.


Introduction
The direct attachment of an artificial limb to the residual stump bone requires that an implant traverse the overlying skin and become firmly osseointegrated (OI) within the medullary canal. This construct provides amputees with an optimum and robust exoprosthetic docking system (Fig 1). The benefits of this technology have been well demonstrated in multiple cohort studies [1][2][3]. Patients with limb loss, who have received percutaneous OI prostheses, have reported vastly improved upper and lower extremity prosthetic function compared to their prior use of conventional socket suspension systems [3][4][5].
Functional ambulation using skeletal docking devices requires that the adjacent endosteal bone in the medullary canal to become strongly osseointegrated within the structure of the implant surface. Over time, this integration depends upon the metabolic process of bone deposition, resorption, and remodeling. These processes occur continually in response to the loading and unloading forces applied to the skeletal bone and are described by "Wolff's Law of bone adaption" [6]. In theory, if OI device design configurations support appropriate loading of the implanted bone, inducing robust bone formation around the implant and avoiding bone loss from "stress-shielding," these percutaneous devices could last the lifetime of the recipient [7,8]. However, some currently used implant systems fail to maintain permanent osseointegration because the design configuration intrinsic to the implant does not support bone deposition and over time encourages stresh shielded bone loss. Such resorption can also result in stomal collapse if the skin around the percutaneous post is attached to this underlying resorbing [9][10][11][12]. To date, there are no techniques that are capable of predicting early implant demise secondary to the failure of the process of osseointegration.
In addition to optimal osseointegration without distal resorption, there are legitimate concerns about the vulnerability of the stoma to superficial and potentially deep infection [2,13]. It has been clinically established that a stable and mature stoma is vital to the longevity of these percutaneous OI devices. Although wound healing processes are known to be impaired around such implants [14], the material surface of the implant in contact with the tissue lining the stoma, must encourage completion of wound healing and avoid a continuous state of chronic inflammation at this interface. Currently, clinical observation and patient-reported outcomes are used to assess the completion of wound healing. As intimated above, objective quantitative techniques that non-disruptively detect the end of bone and stomal wound healing have not been developed. Various invasive techniques have been used in animal models to study the local biological processes involved in these healing tissues [15][16][17]. However, in the clinical setting, due to both ethical and patient safety considerations, serial acquisition of local tissue is difficult, if not, impossible. In contrast, the analysis of whole-blood samples offers a non-invasive alternative that may be used to assess systemic biological processes of interest. Although multiple protein biomarkers have been correlated with bone-related diseases [18], such analyses have not been used as a predictive means for determining early indications of device success, impending failure, or the achievement of the final steady-state of complete wound healing. In contrast to typical assays used for protein detection in tissue samples, RNA sequencing (RNA seq) enables large-scale assessment of tissue activity via measurement of the transcriptome. Thus, RNA-seq was chosen as the preferred method for this investigation. The aim of this study was to use whole-blood RNA-sequencing to identify potential biomarkers involved in inflammation and bone metabolism, which, over time, might indicate the resolution of inflammatory wound healing processes, as well as the progression or completion of osseointegration. Once determined, these mRNA markers may be used to monitor the clinical progression of OI healing using relatively inexpensive PCR techniques.

Study design
Ten veteran unilateral transfemoral amputees, who were recruited to participate in an FDAapproved Early Feasibility Study of a POP device (NCT02720159), were also included in this current study (NCT02564432) with IRB approval (#00073178). The study design, described in detail previously [19], is briefly outlined below. A (Level II) prospective, longitudinal, observational study was conducted where participating patients underwent two separate surgical procedures to fit a transfemoral percutaneous OI implant [DJO global, Austin, TX] with a wait time of 6-8 weeks between procedures. The first procedure involved soft tissue revision, preparation of the residual bone to accept an appropriately sized endo-prosthesis, and device implantation. The second procedure involved coring the skin and muscle covering the implanted device and connecting an exo-prosthesis through the skin opening (i.e., stoma).

Inclusion/Exclusion criteria
Patients were required to be within the US Veteran Health Care System with a unilateral, trauma-related, transfemoral amputation that was not a result of dysvascular disease. Patients were required to have used or be currently using a socket suspension prosthesis and the use of a non-propulsive, passive microprocessor-regulated lower limb prosthesis. Patients had to agree not to participate in high levels of physical activity while in the study.
Exclusion criteria included currently being on active military duty, the loss of more than one, the diagnosis of insulin-dependent or adult-onset diabetes, or a recent tobacco use history.

Sample collection and processing
Fifteen mL venous whole-blood samples were collected from patients at the following time points during the study ( . Single-end 50bp sequencing was performed by the High-Throughput Genomics and Bioinformatics Analysis Shared Resource at the Huntsman Cancer Institute at the University of Utah, which was supported by the National Cancer Institute of the National Institutes of Health (Award Number P30CA042014). Sequence quality was assessed with MULTIQC v1.11 [20] incorporating RNA-seq metrics from FASTQC v0.11.9 [21].

Sequence alignment and feature quantification
Samples were required to have an RNA Integrity Number (RIN) > 6 for inclusion. The resulting sequences were aligned to reference genome GRCh38.94 (STAR [22] via twopassMode).

Count matrix pre-filtering
Because hemoglobin (HGB) sequences can make up a large portion of the RNA contained in whole-blood samples and has previously been shown to mask signals from lower expressed genes [24], HGB counts were identified bioinformatically using Ensembl annotations. Counts for the following HGB genes were removed: HBA1, HBA2, HBB, HBD, HBG1, HBG2, HBE1, HBM, HBQ1, HBZ. In order to be included for further analysis, count filtering was applied, which required at least 2 reads to be assigned to a gene in at least 5 samples.

Differential gene expression analysis
The filtered count matrix was normalized using the default normalization function in DESeq2 v1.32.0 [25] before applying surrogate variable analysis (sva) [26]. The number of surrogate variables generated was five based on the "be" method from the sva package v3.40.0. Pair-wise comparisons were made between PrS1, which served as each patient's baseline, and all subsequent time points using DESeq2 to test for differential expression. For each comparison, surrogate variables, patient, and time points were included in the DESeq2 statistical model. Multiple testing correction was performed using the Benjamini-Hochberg method with a false discovery rate (FDR) of 0.05.  plates were normalized using a common control. Fold-change was calculated using 2 −ΔΔCt values.

Enrichment analysis
The Enrichment analysis and associated visualizations were produced using the ViSEAGO v1.6.0 [27] R package. First, genes were annotated using the gene ontology (GO) [28] biological process annotations. The node size was set to 5, requiring at least 5 genes to be annotated to a GO term, and a Fisher's Exact Test was then applied. Due to the overlap between genes assigned to separate GO terms and that multiple testing methods assume independence of tests, multiple testing correction was not applied in the enrichment analysis. Instead, topGO's elim algorithm v2.44.0 [29] was used. When a node was found to be significant using the elim algorithm, the corresponding genes were marked and removed from all parent nodes prior to applying the Fisher's Exact Test on parent nodes. A p-value below 0.01 was considered significant. To aid in the visualization and identification of relevant GO terms, the semantic similarity between significantly enriched terms was calculated using the method described by Wang et al. [30] Hierarchical clustering was performed and visualized using these similarity measures and the dynamicTreeCut 1.63.0 [31] package. The sequence data has been submitted to the Gene Expression Omnibus (GEO) under accession number GSE194108.

Study design
Ten male adult patients with unilateral trans-femoral amputations were recruited to participate in this IRB-approved study. Two patients experienced device failures. One patient sustained an accidental high-energy injury with a resulting periprosthetic femoral fracture eight months after the Stage 2 surgery. Prior to this unfortunate event, this patient achieved clinically robust osseointegration of the implant and a stable stoma. The second patient was removed from the study early on because of the failure of osseointegration and device loosening. Blood samples were not collected from these patients after the device failure.

Sequence alignment and quantification
One patient was not available for follow-up at PoS2-M9 and after removing samples of low quality, 92 samples were included for downstream analysis (Table 2). On average, the selected samples had a RIN of 7.83. The mean base-pair quality scores were > 32, with an average GC content of 49%. Samples had, on average, 30.7 million reads aligned with 98.7% reads mapped (77.7% uniquely mapped reads, 21% multi-mapped reads). Ribosomal mapped reads were effectively removed during sample processing (0.0058% reads mapped to ribosomal locations).

Gene expression and enrichment analyses
The majority of Differentially Expressed Genes (DEGs) identified (Table 3) were protein-coding genes. Some non-coding genes were identified, which included lncRNA, miRNA, snRNA, and snoRNA gene types. The DEGs highlighted in this section had a log fold-change (log 2 FC) less than -1 or greater than 1 and were grouped into inflammation (Fig 3) and bone-remodeling (Fig 4) categories based on their associated GO annotations. A full list of DEGs identified with their associated time point is included in S1 Table. Over the course of the study, 27 clusters were identified from the enrichment analysis (Figs 5 and 6). Most DEGs were only detected at PoS1 immediately after the first surgery. However, a handful of DEGs was identified at later time points and annotated to inflammation-related GO terms: GRB7, IGHD2-2, IGKV1-16, IL12B, IL1RL2, IL33, NPY5R, PRG2, AND SEMA3 (Fig 3). Two genes, IL1RL2 and IL33, continued to be significant at most time points following the first surgery. Multiple GO terms involved in the inflammatory response and enriched at multiple time points were identified in clusters 13, 15, and 27 ( Fig 5).
Two DEGs involved in bone remodeling were identified: COL2A1 and SOST (Fig 6). SOST was found to be significant at multiple time points. From the enrichment analysis, the GO term, regulation of bone development, was significant at PoS1.

qPCR validation
The qPCR data for fold-change values followed a similar trend to the RNA-seq fold-change values with a correlation coefficient of 0.75 (Fig 7).

Discussion
In this longitudinal study, the systemic expression of RNA from whole blood gene expression was surveyed in a ten-patient cohort, in which all patients were fitted with a percutaneous OI prosthetic system. Repeated measure designs, such as the one employed by this study, decrease the amount of variability experienced between measures by having the patients serve as their own baseline. Overall, the results showed at least four circulating genes of interest that could help to understand early bone remodeling and inflammation: IL1RL2, IL33, COL2A1, and SOST. Multiple DEGs, annotated to inflammation-related GO terms, were identified from the enrichment analysis. Most DEGs were only significant at PoS1. However, two genes, IL33 and IL1RL2 were consistently up-regulated following the first procedure. Interestingly, IGHD2-2 was significantly down-regulated at time points when IL33 returned to baseline (PoS2-M3 and PoS2-M9), the reason for this finding is not readily apparent, but is currently being investigated. IL-33 and IL1RL2 (the latter which serves as a receptor for IL36) are actively involved in the inflammatory process. In mice, IL-33 has been shown to be involved in bacterial defense by inducing nitric oxide production [32]. In another study, inflammation was significantly reduced in IL-36 -/mice who were exposed to Staphylococcus aureus (S. aureus) [33]. Interestingly, our previously reported study investigating the microbiome of the stoma in these same patients identified an increase in the relative proportions of Streptococcus and Staphylococcus bacteria over time [19]. This result indicates that the periprosthetic tissues may be increasingly exposed to colonization by S. aureus, and that the elevated levels of IL33 and IL1RL2 may contribute to providing local immunity against superficial infection. Conversely, in the absence of a matured stoma and a skin seal against the implant, it is possible that the sustained upregulation of IL-33 and IL1RL2 may indicate increased ingress of potentially pathogenic but colonizing bacteria into the peri-prosthetic tissue, albeit in the absence of clinical infection. It is also possible that increased systemic expression of IL-33 and IL1RL2 could be attributed to the patient's adaptive immune responses to minimize infection associated with the infiltration of bacteria into the periprosthetic tissue from the percutaneous stoma. However, further investigation is warranted to confirm these assumptions.
COL2A1 is typically expressed in cartilaginous tissue and is involved in the production of collagen type 2 extracellular matrix [34]. For this gene, we observed a down-regulation of over 20-fold in the circulating blood at PoS2-M3. Immediately following the second surgery, our patients were able to load their endoprosthesis. Post-loading, we expect the host bone to undergo accelerated remodeling as supported by our ovine preclinical data, which showed an accelerated bone remodeling period within 3-months of post-implantation [35,36]. Furthermore, COL2A1 expression has been reported within the local bone tissue during the remodeling process [37]. The question remains as to whether or not the downregulated Col2A1 at 3-months following the second surgery marks the end of the accelerated bone remodeling period. Although COL2A1 is not directly related to the bone matrix, scientists increasingly recognize the emerging role of chondrocyte-to-osteoblast transdifferentiation in the endochondral ossification process [38]. Moreover, in a mouse dental implant study, COL2A1 was upregulated during the osseointegration process within the local oral tissue [39]. Perhaps, this transdifferentiation is also the molecular pathway plays a role in the osseointegration of OI implants. This concept also warrants further investigation.
In the presented dataset, SOST expression levels were significant and cyclic following the first surgery. This may indicate control of local bone resorption/deposition and possibly ongoing changes in the rate of osseointegration. Thus, SOST is another potential gene of interest. Although SOST mRNA is secreted primarily by osteocytes, very small quantities are reported to be secreted by basophils, eosinophils, and neutrophils [40]. The protein encoded by the SOST gene, sclerostin, belongs to the family of bone morphogenetic protein (BMP) antagonists and is known to regulate bone turnover, i.e., negatively regulate bone formation [41]. Sclerostin is also involved in regulating osteocytes, which are major mechano-sensors in the bone. In our study, SOST expression was initially upregulated between the first and second surgeries. This may indicate bone resorption due to nonloading during this period which is when patients are advised to not load their implant limb. This phase was followed by decreased SOST expression, with loading between second surgery and PoS2-M3. A similar result was observed, in a mouse study, within the local tissue [39]. Our previous in vivo ovine studies supported this observation in which the growth of new bone into the cancellous titanium structure at the endosteal interface occurred very rapidly but required over 12-months to reach a steady state [35,36,42]. It is also known that differentiating osteoblasts at the end of the bone deposition period are known to transition into osteocytes. Thus, at the end of the accelerated bone remodeling period, this transitioning phase may have resulted in the observed SOST overexpression. Furthermore, multiple studies in the literature highlight the importance of sclerostin in bone remodeling and development. In a SOST knockout mouse study, Shu and his coworkers have shown that inhibition of this gene improved implant fixation in osteoporotic bone [43]. Other authors have shown that when the osseointegrated implant was loaded, local SOST expression was downregulated [44]. Thus, SOST serves as another serological gene candidate for studying osseointegration and bone-resorption following any orthopedic implantation surgery. The qPCR fold-change values of these four genes, in a subset of patients, followed a similar trend when compared to the RNA-seq log-fold-change values with a correlation of 0.75 validating our DEG analyses data.
The current study is not without limitations. Although the samples had, on average,~30 million reads, this depth of sequencing coverage may be inadequate with the removal of HGBrelated counts and the presence of multiple different cell types contained in the blood. To alleviate these issues in future studies, one may incorporate a globin mRNA depletion step to allow for wider read coverage and/or increase the depth of coverage per sample. Another limitation included the small number of patients enrolled in the study. Ten patients were initially enrolled, including two of whom which experienced implant failures. Additionally, it is possible that the statistical model used for identifying DEGs may suffer from overfitting by including multiple surrogate variables in the model. The majority of DEGs were detected at PoS1, and most genes approached their initial baseline expression shortly following the second surgery. However, a handful of inflammatory-related genes continued to be significant at twelve months following the Stage 2 Surgery. It is impossible, with the current study design, to determine if these genes would return to baseline or continue to be up-regulated given additional time. Although four genes were validated using qPCR, further protein expression validation is warranted. Such assays may involve collecting local stomal exudate (which can, but not always, be present clinically) or local tissues where respective proteins may have been translated rather than in the blood. Future animal studies and human trials will be used to confirm these findings.
In conclusion, this preliminary study of whole-blood RNA-seq from 10 transfemoral amputees, who had previously undergone percutaneous osseointegrated skeletal docking of a prosthetic lower limb, uncovered two major results 1) Two inflammatory-related genes (IL33 and IL1RL2) were found to be differentially expressed by 20-fold. These genes may serve as important biomarkers in the serial assessment of stomal wound healing. 2) Two other biomarkers related to bone remodeling (COL2A1 and SOST) were identified. If verified by future studies, these four genes may serve as markers for predicting optimal bone remodeling and stomal tissue healing following the percutaneous OI implantation procedure.
Supporting information S1